order of included R functions:

• lm()
  
• summary()
  
• ggplot() + statqq()
  
• ggplot() + geom_smooth()
  
• tidy(), augment(), glance()
  
• summ(), plot_summs()
  
• autoplot()
  
• glm() (different families)
  
• lmer()

• Linear Regression
• Multiple Linear Regression
• Generalized Linear Models
• Briefly: Mixed effects model

In a linear regression, we aim to find a model:

• that represents our data and
• can give information about the association between our variables of interest. The command in R for a linear model is
  

lm(y ~ x, data = data).

This is called the formula notation in R.

• The formula-type code with the data assignment at the end might seem non-tidy.
• But, we can use the pipe %>% in the following way: data %>%
lm(y ~ x, data = .)

• with the “.” we assign the data in the pipe

Let’s first have a look at the summary table of the example data set, by using the summary() command:

##       day1             day2               day3      
##  Min.   :-1.611   Min.   :  0.2396   Min.   : 0.00  
##  1st Qu.: 3.031   1st Qu.:  9.9763   1st Qu.:13.25  
##  Median : 5.619   Median : 15.9973   Median :19.37  
##  Mean   : 5.428   Mean   : 16.7252   Mean   :20.27  
##  3rd Qu.: 7.772   3rd Qu.: 20.6984   3rd Qu.:28.05  
##  Max.   :12.142   Max.   :143.8187   Max.   :39.03  
##                   NA's   :4          NA's   :16

Before we start with the linear regression model, we need to get an idea of the underlying data and its distribution. We know that the linear regression has the assumtptions:

• Linear relationship
• Multivariate normality
• No or little multicollinearity
• No auto-correlation
• Homoscedasticity

In the plots, we could already see, that day 2 and day 3 seem to be associated. This is obvious when drawing a line in the plot. Let’s now perform a linear regression model in R.

lm(day2 ~ day3, data = tki_demo)

• As said before, the first argument in the code is y, our outcome variable or dependent variable. In this case it is day2.

• The second Argument is x, the independent variable. In our case: day3.

• We also specify the data set that holds the variables we specified as x and y.

Now we want to look at the results of the linear regression. So how do we get the p-value and \(\beta\)-coefficient for the association?

In R, we add the summary() function to the lm() function, like so:

summary(lm(y ~ x, data = data))

We can also store our model results in a variable:
model1 <- lm(y ~ x, data = data), and then use summary: summary(model1)

• The output before contains a lot of relevant information, but it is not straighforward to access the individual parameters like p-values and betas.
• The broom R package is in line with the “tidy” data handling in R and turns the linear model results into an easy accessible tibble format:

tidy(lm1)

## # A tibble: 2 x 5
##   term        estimate std.error statistic    p.value
##   <chr>          <dbl>     <dbl>     <dbl>      <dbl>
## 1 (Intercept)   4.52      0.887      5.10  0.00000214
## 2 day3          0.0382    0.0402     0.951 0.345

• Often we want to visualise the coefficients in the model to see their impact on the outcome, or visualise the coefficient of specific variable in two models, that differ only in the adjusted covariates.
• The jtools package has a nice function to do this very easily, utilising ggplot2:
  plot_summs()

plot_summs(lm1)

plot_summs(lm1, lm2)

plot_summs(lm1, lm2, coefs = "day2")

• The generalized linear model works very similar to the lm() function. We use the glm() function instead.
• The lm() function always assumes a linear distribution, whereas in the glm() function, we have to assign the distribution (“link”) and a family type ourselves.
  
• For a logistic regression, we have to assign: family = binomial(link = “logit”)
• For a poisson regression, we assign: family = poisson(link = “log”)

When we want to also account for random effects or clusters such as intervention if the effect of this is not what we are interested in, or study center in the case of a multi-center study, we use the lme4 package.

The formula works in the following way:

• lmer(y ~ x + (1 + x | randomEffect), data = data) , for random slope
  
• lmer(y ~ x + (1 + x | randomEffect) + (1 | otherVariable), data = data), including otherVariable as a variable that has an impact on the slope

library(lme4)
lmer(day1 ~ day2 + (1 + day2 | intervention), data = tki_demo )

## Linear mixed model fit by REML ['lmerMod']
## Formula: day1 ~ day2 + (1 + day2 | intervention)
##    Data: tki_demo
## REML criterion at convergence: 496.8176
## Random effects:
##  Groups       Name        Std.Dev.  Corr 
##  intervention (Intercept) 1.278e-04      
##               day2        5.528e-06 -1.00
##  Residual                 3.147e+00      
## Number of obs: 96, groups:  intervention, 3
## Fixed Effects:
## (Intercept)         day2  
##     4.79533      0.03436  
## convergence code 0; 1 optimizer warnings; 0 lme4 warnings

• t-test:
  • t.test(datx, daty)
• wilcoxon test:
  • wilcox.test(datx, daty)
• chi-square-test:
  • chisq.test(group1, group2)

• Linear regression and multiple linear regression both work by using lm()
• Models with for other distributions (e.g. logistic regression):
  • glm()
  • specify family and
  • link
• Linear mixed effects model: lmer()

• All models work with the broom and jtools packages
• tidy results: tidy()
• more results: e.g. augment(), glance()

• With the jtools package
• export_summs() for regression table
• plot_summs() for forest plots